en 






< 



o 



a 



AN IMPROVED TALBOT METHOD FOR NUMERICAL LAPLACE 
TRANSFORM INVERSION 

BENEDICT DINGFELDER AND J.A.C. WEIDEMAN 



^-^ Abstract. The classical Talbot method for the computation of the inverse Laplace 

^vj transform is improved for the case where the transform is analytic in the complex 

plane except for the negative real axis. First, by using a truncated Talbot contour 
rather than the classical contour that goes to infinity in the left half-plane, faster 
convergence is achieved. Second, a control mechanism for improving numerical 
stability is introduced. These two features are incorporated into a software code, 

j»^ whose performance is assessed on transforms from tables as well as from actual 

applications. It is shown that even when the transform has singularities off the 

I I negative real axis, rapid convergence can still be achieved in many cases. 

B 



1. Introduction 



Let f{t) be real-valued and piecewise continuous on [0, oo) and of exponential 
order / = 0{e'^°*), t -^ oo, for some real constant 70. Then the Laplace transform 
of /(f) is defined by 

^"^ /"OO 

> f (z) = / e-^*f(t)dt, Rez>7o. 

I/-) Jo 

^^ This paper deals with the numerical solution of the inverse problem, i.e., given 

_J F{z), compute f{t) at a specified value of t. The function F{z) is analytic in the 

half -plane Re z > 70, and if it can be evaluated in the complex plane (as opposed 
to just on the real axis), this analyticity can and should be exploited by good 

7^ numerical methods. One such method is numerical quadrature applied to the 

7—i inverse formula 

. ^ /(O = ^ r^'^ e^'f (z) dz, Re 7 > 70, (1) 

X 



Z7TI J'y — ioo 



;_i known as the Bromwich integral. Two effective quadrature rules are the simple 



trapezoidal and midpoint rules, used in combination with contour deformation. 
The purpose of the contour deformation is to exploit the exponential factor 
in (1). In particular, assume the path of integration in (1) can be deformed to a 
Hankel contour, i.e., a contour whose real part begins at negative infinity in the 
third quadrant, then winds around all singularities and terminates with the real 
part again going to negative infinity in the second quadrant. Examples are shown 
in Figure 1 below. On such contours the exponential factor causes a rapid decay, 
which makes the integral particularly suitable for approximation by the trapezoidal 
or midpoint rules. The contour deformation can be justified by Cauchy's theorem, 
provided the contour remains in the domain of analyticity of F(z). Some mild 
restrictions on the decay of F(z) in the left half-plane are also required; for sufficient 
conditions, see (Talbot 1979, Trefethen, Weideman and Schmelzer 2006). 



2 BENEDICT DINGEELDER AND J.A.C. WEIDEMAN 

Suppose such a Hankel contour can be parameterized by 

C: Z = z{e), -TT^ ^ 7T, (2) 

where Rez(±7r) — —00. Then 

We approximate the latter integral by the N-panel midpoint rule with uniform 
spacing h = 2ti/N, which yields 

1 N 1 

/(O « ^ E'^'^'^^'f(^(^fc))z'(efc), 0fc = -n+{k--)h. (4) 



If the contour (2) is symmetric with respect to the real axis, and if F{z) = F{z) 
(which is the case for /(f) real-valued), then half of the transform evaluations can be 
saved. That is, one needs to consider only quadrature nodes in the upper (or lower) 
half-plane. When comparing numerical results it is important to keep this point in 
mind as several other papers, including (Duffy 1993, Weideman 2006, Weideman 
and Trefethen 2007), considered quadrature rules with a total of 2N nodes. 

Popular contours C are the parabola (Butcher 1957, Gavrilyuk and Makarov 
2005) and the hyperbola (Sheen, Sloan and Thomee 2003, Lopez-Fernandez and 
Palencia 2004, Gavrilyuk and Makarov 2005), as well as the cotangent contour 
introduced in (Talbot 1979). Here we consider the Talbot contour in the form 

N 
z{e) = —(,{0), ^{e) = -cr + iiecot{ae)+vie, -Tr^e^/r, (5) 

where a, ji, v and a. are constants to be specified by the user. In the original 
Talbot contour the parameter a did not appear, i.e., a = 1. For reasons outlined 
below, the modification we propose is to consider < a < 1. The scaling factor 
N/t was introduced in the original paper (Talbot 1979) and also proposed in 
(Weideman 2006). Assuming a fixed value of t, this is the appropriate scaling for 
maximizing the convergence rate in the trapezoidal /midpoint approximation as 
N ^- 00. 

In the original Talbot contour (a = 1), it follows that Rez(±7r) = —00. If the 
trapezoidal rule is used, this means the point at infinity is a quadrature node 
and the corresponding function value is set to zero in the sum (4). As this seems 
wasteful, we recommend the use of the midpoint rule instead. 

The practical issue addressed in this paper is the choice of the parameters 
(J, n, V and a in (5). The main considerations are summarized in Figure 1. The best 
contour has to strike a balance between passing too close to the singularities, or 
too far from them in which case the exponential factor in (3) becomes too large. 
The quadrature nodes on the best contour should also extend just far enough into 
the left half -plane to reach the desired accuracy but their contributions should not 
be negligibly small. 

In the case a = 1, suggestions for parameter values w^ere made in the original 
paper (Talbot 1979). These were implemented in ACM TOMS Algorithm 682 (Murli 
and Rizzardi 1990). The performance of this algorithm on a set of challenging test 
problems was assessed in (Duffy 1993). In Section 4 we show that the modified 
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Figure i. Three possible Talbot contours (5), each shown with N = 24 
nodes in the midpoint rule. The middle contour is close to optimal, in 
the case where t = \ and the singularities of ¥{z) are located on the 
negative real axis. Its outermost conjugate pair of nodes reaches right 
up to the dash-dot curve, where the magnitude of the exponential factor 
in (3) becomes less than the desired accuracy (here taken to be machine 
precision, e w 2.2 x 10^^^). The inner contour is suboptimal because (a) 
its outermost two nodes do not reach quite up to the dash-dot curve and 
(b) it passes too close to the singularities. The outer contour is suboptimal 
because (a) the contributions of the outermost pair of nodes are less 
than the desired accuracy and therefore wasteful, and (b) the exponential 
factor in (3) is too large in the right half-plane. 



contour (5) with the parameter choices proposed here can produce the same 
accuracy but often with significantly fewer terms in the sum (4). 

It is unrealistic to try and find optimal contours for all possible transforms 
F(z). That is why we limit the analysis here to the situation that seems to occur 
most frequently in applications, namely the case where all singularities of F{z) are 
located on the negative real axis. We shall see, however, that the contour parameters 
derived here also work reasonably well when singularities are off the real axis (but 
not too far, depending on the value of t). 

In the case a = 1, an optimized set of parameters a, pi, v was proposed in 
(Weideman 2006). When all singularities of F(z) are located on the negative real 
axis, the achievable convergence rate is 0{e^^^'^'^'^^) as N ^> 00 for fixed f > 0. 
Because Re2(±7r) = —00 in the original Talbot method, however, too many nodes 
end up far out in the left half-plane where they make no contribution to the 
midpoint sum (4). We acknowledge (Trefethen 2006) for this observation, and the 
suggestion to consider a truncated Talbot contour instead, with < a < 1. In this 
case we refer to (5) as the modified Talbot contour. 
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Optimal parameters a, j,i, v and tx. for the modified contour were derived by one 
of the authors of the present paper, and reported in (Trefethen et al. 2006), later to 
be included in text books such as (Strang 2007, Sect. 5.3). Details of the derivation 
have not been published, however. Having simplified and improved the derivation 
recently, we now present the details here. 

The freedom that the parameter a offers means that the convergence rate 
0(-g-O.949N) fQj. t^g g^gg fl, = 1 can be improved to O [e-'^-^^^^ ) , N ^ 00, t fixed. 
This is assuming singularities of F(z) on the negative real axis. By comparison, with 
the same assumption on the singularity location, the convergence rates achievable 
by the parabola and hyperbola are, respectively, (!)(e~^'^^'^^) and 0{e~^'^'^^^); see 
(Weideman and Trefethen 2007). 

It should be noted that the scaling factor N /t in (5) has a few negative conse- 
quences. First, because the nodes z{6i^) depend on N, function values cannot be 
re-used as N is increased. The increased convergence rate compensates for this, 
however. Second, for fixed t the contour can move far into the right half-plane as N 
is increased. Because of the exponential factor in (1), the terms in the summation 
can become large and hence suffer from floating-point cancelation errors. We 
address this issue in Section 3. Third, because the quadrature nodes depend on t 
as well, the same contour cannot be used for two or more distinct values of f . In 
this case the method would be inefficient when the transform F{z) is expensive 
to evaluate. We shall not address the latter issue here, but note that for the case 
a = 1 it was considered in (Rizzardi 1995) and for the parabolic and hyperbolic 
contours it w^as considered in (Weideman and Trefethen 2007). 

The outline of the paper is as follows: In Section 2 we present the main error 
analysis that leads to the proposed values of a, pi, v and a. The analysis is similar 
to the one given in (Weideman 2006). It is slightly more intricate because of the 
additional parameter a, but at the same time w^e have made some improvements 
that simplify the analysis. A strategy for improving numerical stability for large 
N is discussed in Section 3. For now it only makes provision for roundoff error 
control when the singularities of F(z) are located on the negative real axis. Section 4 
contains the numerical tests, where we compare the results of the modified contour 
(5) with parameters proposed here, with Algorithm 682 (Murli and Rizzardi 1990). 
Further tests and the details of a code, that implements the proposed method are 
discussed in Section 5. The code, written in MATLAB (MATLAB 2012), is publicly 
available for download from the web site (Dingfelder and Weideman 2013). 



2. Error estimates 

Our basic error estimate is a well-known contour integral formula for the error 
in the trapezoidal rule approximation of functions that are analytic, periodic, and 
real-valued; see for example (Kress 1998, Thm. 9.28). We retain the property of 
analyticity but extend the theorem as follows: (A) we consider non-periodicity, 
(B) we consider complex-valued functions, and (C) we consider the midpoint rule 
rather than the trapezoidal rule. (A) is necessary because the integrand in (3) is 
not periodic. (B) is necessary because the integrand is complex-valued, and has 
very different analytic properties in the upper and lower half-planes. (C) is not 
really necessary, only a matter of convenience. All three of these modifications can 
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be readily incorporated into the proof given in (Kress 1998, Thm. 9.28), which is 
an elementary application of the residue theorem. We omit the details. 

Theorem 2.1. Let d^^ be defined as in (4). If g : D ^ C is analytic in the strip D = {0 e 
C:— 7r<Re0<7r and — d <lra 6 < c} where c,d > 0, then 

N 



J-n l>i ,,_i 



Here 



and 



E+(7) = ^ / +/ +/ ){i + itM^))g{0)de 

•^ \J — 7r J — TT+i'y Jn+ij/ ^ 

1 / r-n-iS rn-iS rn \ \jfl 

for all < J < c, < 5 < d and N even. For odd N, replace the tan(N0/2) ivith 
-cot{Ne/2). 

The contours of integration for computing the error integrals E+(7) and E_((5) 
are shown in the left diagram of Figure 2. We remark that if g{6) is real-valued 
on the real axis, then g{9) = g{6) and d and c can be taken to be equal. When it 
is moreover 27r-periodic, the contributions on the sides Refl = ±7r cancel. In that 
case, if we define E(7) = £+ (7) + E_ (7), then 

E(7)=Re/ {l + itan( ))g{e)de. 

J — TT+i'y ^ 

Here, and below, we assume N even but the analysis for odd N is similar. By 
analyzing the behavior of the tangent function in the complex plane, this can be 
bounded by (Kress 1998, Thm. 9.28) 

,^, ,, 47rM 

where M is a bound on the magnitude of g in the strip D. For a given value of N, 
two quantities control the size of the error: c, which is the width of the strip of 
analyticity, and M, w^hich is governed by the growth of g in the complex plane. 
Here w^e apply the error formula of Theorem 2.1 to the case 

g{e) = ^.e-^'^'V(z[Q))7!{Q), -n^e^n, (6) 

where F(z) is the transform to be inverted. This is not a real-valued function, 
which is why the rectangle of integration in Figure 2 is not shown symmetric 
with respect to the real axis. As we shall see (and similar observations were made 
in (Weideman 2006, Weideman and Trefethen 2007)), in the upper half-plane the 
width of the strip is restricted by the singularities oi F(z{6)), while in the lower 
half-plane the restriction is the size of e^(^''. 

To keep the analysis tractable, we limit our discussion primarily to the model 
transform 

F{z) = {z + \)-\ A>0, (7) 

i.e., singularities on the negative real axis only. In Section 4 we shall consider a 
matrix version of this transform when considering applications to PDEs. 
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Figure 2. Left: Rectangular contour in the complex 6-plane for evaluating 
the error integrals E+ (7) and E_ (S) of Thm. 2.1. Right: Deformed contour, 
along which the magnitudes of the integrands of E+ (7) and E_ (S) are 
approximately constant. The dash-dot curve in both diagrams shows a 
possible image of the negative real axis of the z-plane, where singularities 
might be located. 



For (7) and similar transforms, we can estimate the error using Theorem 2.1 
and the rectangular contour shown on the left in Figure 2. For a better estimate, 
one can try and deform the contour as shown on the right in Figure 2. 

To explain how this can be done, consider the error integrals in Theorem 2.1. 
Ignoring the constant factors 1/2 in those integrals, as well as the factor 1/ {Irci) 
in (6), the error formulas can be expressed as 



E± 



h±{e) 



.Nh±{e) , 



c± 

m 



N- 



,Ne 



log (1 ± f tan( — )) + logF(2(0)) +logz'(0) 



where C+ and C_ refer to contours in the upper and lower half-planes, respectively. 
Our strategy is now as follows: suppose it is possible to deform the rectangle on 
the left in Figure 2 to a contour along which Reh±{6) = —c for some positive 
constant c. Then 



C± : Reh±[e) 



|E±|=0(e 



-cN\ 



(8) 



and the value of c can be maximized by picking the best among all such contours 
C±. 

The function h±{9) defined above is too complicated to analyze in this manner, 
and we make a few approximations in the case N ^ 1. First, we drop the terms 
involving F{z{Q)) and z' [9). This means F(z(0)) should not be zero nor exponen- 
tially small on any part of the contours C±, and likewise z' (Q) 7^ 0. Regarding the 
term involving the tangent function mh±{9), we shall approximate it for N 2> 1 
as follows. By converting the tangent function to exponentials one obtains, with 
6 = X + iy, 

l±i\an[ — {x + iy)) = 



Ny^^iNx, 



■62 



%. 



ilN.T 
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Taking the logarithm yields 

MO 
log (1 ± f tan( — )) ~ ^Ny, N ^ oo, (9) 

valid for fixed, nonzero y. The upper and lower sign choices correspond to y > 
and y < 0, respectively. Accordingly, we approximate the contour Reh±{6) = —c 
in (8) by 

C± : Re l,{x + iy) =F y = —c- (10) 

First, consider the upper half of the 0-plane. As mentioned above, the error 
contour C+ is restricted by the location of the singularities of F(z(0)) in the upper 
half of the 6-plane. For now, consider the model problem (7). Its pole is located 
at z = —A, or f = —At/N. In keeping with the limit N ^- 00, f fixed, we consider 
therefore ^ = 0. The same argument will apply if there is a branch-cut on the 
negative real z-axis, for in practice we need to consider only a finite section of the 
branch-cut (e.g., the section between the origin and the dashed line in Figure 1). 

We therefore need to examine the zeros of (,{0), i.e., 

-cr + }i 9 cot{a9) +vie = 0. 

It is no easy task to examine the complex roots of this equation for all possible 
parameter values, so we resorted to numerical experimentation, restricting searches 
to the strip — tt ^ Re ^ n, which is of interest here. For certain parameter choices 
there are roots in the lower half of the 0-plane. For other parameter choices 
the roots are in the upper half-plane, located symmetrically with respect to the 
imaginary 0-axis. For yet other parameter values the roots are located precisely 
on the imaginary 0-axis, with the one closest to the real axis the critical one that 
limits the contour C+. 

Another limiting case is the fact that ^'(0) may not be 0, as remarked below (8). 
By numerical experiments, we came to the same conclusion as (Weideman 2006), 
namely the configuration that gives the widest domain of analyticity in the upper 
half-plane occurs when these two critical points coincide precisely on the positive 
imaginary 0-axis. (This is the point marked P in Figure 2.) This leads to the two 
equations l,{iy) — 0, ^'(fy) = for some y > 0. By considering (10), we conclude 
that y = c. If we further let ^> ± tt on the contour C+, we get ^(±7r) = — c, which 
is by symmetry just one equation, not two. In summary, in the upper half-plane 
we require that 

C+: C(fc)=0, r(;c)=0, l[n) = -c. (11) 

By using these three equations, it is possible to solve for (y, }i, v) in terms of 
[oi,c), namely 

(T = loL(p-B, /( = 2sinh^(ac)B, i/ = (sinh(2ac) - 2ac)B, (12) 

where 

csrn^(a7r) 

^ T 7 2 ■ 

2.cc(P-svc\ {oiTi) — TT sin(2a7r)sinh (ac) 

It remains to fix the parameters a and c. For this, we consider the lower half 

of the 0-plane. Here the contour C_ is not restricted by analyticity but by the 

size of the exponential factor in (3). Again numerical experimentation confirmed 

the findings of (Weideman 2006), namely that there are two saddle points in the 
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Figure 3. Three possible deformations of the rectangle of Figure 2. The 
middle one maximizes c, ^vhich is both the value of the intersection on 
the positive imaginary axis and the decay rate in the error estimate 

£ = C(e-';N). 



lower half of the 0-plane that are critical. (These are the points marked S and S' 
in Figure 2.) By making C_ pass through these saddle points, we can enclose a 
region as big as possible. 

Denote the saddle point locations by = Xg + i\js- By inserting this into (10), 
and taking partial derivatives, yields 

C_ : Re ^(xs + ixjs) + i/s = -c. Re ^'(xs + i\js) = 0, Im ^'(x,- + iys) = 1- (13) 

We insert the expressions for a, \i and v from (12) and fix a value of ot. Assuming 
a solution exists, the three equations (13) can then be solved numerically for the 
unknowns x^, j/s, and c. By varying a, the value of c can thus be maximized. 

We made no attempt to establish theoretical existence of a solution of the 
equations (13), but direct numerical computation indicates that such a solution 
exists for a in the interval [0.51,0.82], roughly. Some of these configurations are 
shown in Figure 3, with the middle one corresponding to the maximum value of c. 
We computed this value of c using a univariate optimization code, and found 



0.6407; 



1.3580. 



(14) 



Substitution of these two values into (12) yields the modified Talbot contour 

z(0) = —(-0.6122 + 0.5017 0cot(O.64O70)+ 0.2645/0), -n<^Q<^n, (15) 

as reported in (Trefethen et al. 2006). For completeness, we record that the corre- 
sponding saddle points are located at Xg + iy^ = ±3.4208 — 2.3438 i, correct to all 
digits shown. 

An alternative contour in which the cotangent is replaced by a rational function 
was respectively suggested and analyzed in (Talbot 1979) and (Weideman 2006). 
Applying the same strategy as outlined above, we derived 

3.0232 02 



N 



z(0) = -(0.1446+ ^,-3^^j^^ 



- 0.2339 f( 



-71^ ^ TT. 



(16) 



This rational contour is virtually as good as the cotangent contour, as it gives a 
convergence rate of 0{e~^-^^^^), compared to the 0{e~^-^^^^) of (15). With either 
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contour, an accuracy of close to machine precision e ~ 2.2 x 10^ can be reached 
for values of N around 26 or 28, i.e., about 13 to 14 transform evaluations. 

We emphasize that we have assumed here that all singularities of F{z) are 
located on the negative real axis, and f (z) is not exponentially small as discussed 
below (8). In Section 4 below, we test how well the contour (15) works when these 
assumptions are violated. But first, we address the issue of roundoff error that can 
contaminate the solution when N is too large. 



3. Roundoff error control 

When N gets large, the contours (15) and (16) move too far into the right half- 
plane. The exponential factor in (4) gets too big, and there is loss of precision in 
the summation. Using the same arguments as presented in (Weideman 2010) for 
the parabolic contour, we model this error by 

K = C'(ee^(0)') =O(ee^?(0)), (17) 

w^here e is the unit roundoff. 

Let the implied constants in (8) and (17) be denoted by ki and k2, respectively. 
Then the roundoff error becomes significant when 

k^e'"^ = k2e 6^^'^°^ => c + ^(0) +N-i log(e/fco) = 0, (18) 

where we defined fcg = k-[/k2- 

To check the accuracy of (18), consider e ~ 2.2 x 10~^^ and assume for the 
moment that fcg w 1. With the values of the contour parameters as defined by (12) 
and (14), a numerical solution of (18) yields N ~ 24 (25 in the case of the rational 
contour (16)). An examination of the convergence curves shown in Section 4 will 
confirm that N = 24 is often close to the critical value where roundoff errors start 
to dominate. 

When k() is not 0(1), this critical value of N may be larger, say N = N*. The 
value of N* can be estimated by analyzing the convergence behavior. For example, 
we assume that for N ^ N* the error behaves like f{t) — fj^{t) = kie^'^^ , where 
/js/(t) denotes the approximation (4). By computing /n(0 for a range of values 
of N, it is straightforward to detect at what value of N the convergence turns 
from exponential decay to exponential growth. By inserting this value of N* into 
(18), and using the values of the contour parameters defined by (12) and (14), it is 
possible to compute an approximate value of fcg. 

We use this value of ko to compute new contour parameters for values of 
N > N*. Note that ^(0) = —(t+ }i/a, which depends on c and oc via (12). This 
means the equation in (18) has two unknowns, so we fix oc at the value (14). Hence 
c can be solved by applying a univaritate rootfindrng routine to equation on the 
right in (18), and following that the new contour parameters are computed from 
(12). 

We can justify this approach by looking at the asymptotic behavior of the 
parameters. Letting N ^ 00 in (18) yields 

c = - log(e/fco)N-i + 0(N-3), ^(0) = 0{N-^). 



10 BENEDICT DINGEELDER AND J.A.C. WEIDEMAN 

Substitution of these two expressions into the error estimates (8) and (17), respec- 
tively, shows that both reduce to 0{e). Hence errors can be expected to remain at 
the roundoff level for large N. 

Because ^(0) = 0{N^^), the contour is not only prevented from moving too far 
into the right half-plane, it actually moves back to the left (and becomes narrower) 
when N > N* . For this reason this roundoff control strategy only works when the 
singularities of F{z) are strictly on the negative real axis. In the next section, w^e 
investigate how w^ell it works in practice. 

We remark that roundoff control for the parabola and hyperbola was proposed in 
(Weideman 2010) and (Lopez-Fernandez, Palencia and Schadle 2006), respectively. 
An entirely different strategy was developed by one of the authors of the present 
paper, based on methods that automatically compute a contour of optimal stability 
for (1). Details can be found in (Dingfelder 2012) and a forthcoming publication 
on this topic. 

4. Applications and Comparisons 

In this section we verify the theoretical results of the previous sections by 
applying the method to various model problems. These include a matrix transform 
that can be used for the time integration of parabolic PDEs, as well as some scalar 
transforms. The latter include four of the transforms considered in (Duffy 1993), 
all of which were taken from actual applications. In all these cases the inverses 
of the transforms are explicitly known as integrals or series, which enabled us to 
compute errors. Some of these transforms have singularities exclusively located on 
the negative real axis. The other problems deal with singularities on the imaginary 
axis and on the positive real axis as well as with various branch cuts in the complex 
plane. We shall see that good convergence is still possible for these cases, despite 
the fact that the parameter analysis of Section 2 assumed singularities strictly on 
the negative real axis. 

In all figures, the relative error is plotted against the number of nodes, N, 
in the midpoint approximation (4). For reference, we also plot the theoretically 
predicted convergence rate 0(e~^-^^^^), as a dash-dot curve. We also remind the 
reader about the N vs. 2N convention mentioned below (4). When compared to 
the error curves in (Duffy 1993, Weideman and Trefethen 2007, Weideman 2010), 
for example, the values of N on the horizontal axes of the graphs below should 
therefore be halved. 

4.1. Singularities on the negative real axis. We start with the model problem (7) 
and its matrix analogue 

Fi(z;A) = (z + A)"\ Ft{z;A) = {zI + A)-\ (19) 

Here A > and A is a symmetric positive definite matrix, with I the identity 
matrix of the same size. 

The matrix version arises, for example, when semi-discretizations of parabolic 
PDEs are solved by Laplace transform techniques (Sheen et al. 2003, Lopez- 
Fernandez and Palencia 2004, Gavrilyuk and Makarov 2005). Consider 

Uf + Au = =^ 2U(z) -uo + AU(2) = => U(z) = Fi(z;A)uo, 
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where A is an M x M matrix representation of the negative Laplacian, u(t) an 
M X 1 vector of unknowns with U(z) its Laplace transform, and uq the initial 
condition. The inverse formula (i) yields 

u(f) =exp(-AOuo = —Je^'{Fi{z;A)-ao)dz. (20) 

In this case the singularities of the transform are the negatives of the eigenvalues 
A of A, which are located on the negative real axis because of the assumptions 
on A. The right-hand side is approximated by the midpoint sum (4), where each 
quadrature node requires the solution of a linear system with coefficient matrix 
[z{9if)I + A) and right-hand side uq. 

Our test example is the heat equation Ut — 0.01 V^m = on [0, 1] x [0, 1], sup- 
plemented with homogeneous Dirichlet boundary conditions. We take A as the 
familiar block-tridiagonal matrix based on the 5-point finite difference approxi- 
mation to the negative of the Laplacian, which is known to be positive definite. 
The right-hand side uq is taken to be random, and the reference solution u(f) was 
computed using the spectral decomposition of A, which is explicitly known. 

The results for this example are shown in Figure 4. Note that we have plotted 
the error here in the computed value of exp(— Af) uq, i.e., only the temporal error 
The error in the actual solution of the PDE is not shown but it would not reach 
such small values unless a super accurate space discretization is used. 

Figure 4 confirms that (a) the modified Talbot contour converges as predicted 
by the error analysis of Section 2, namely 0{e^^-^^^^), and (b) the roundoff control 
strategy described in the previous section works well in practice. Ten digit accuracy 
is reached with N = 16 or fewer nodes, i.e., no more than 8 transform evaluations 
are required (each of which involves the solution of a sparse linear system). Almost 
full accuracy is reached for all N > 24. 

The remaining transforms in this section are scalar problems from actual appli- 
cations, as collected in (Duffy 1993). As an example of a transform with a series of 
poles on the negative real axis we consider 

^ (lOOz-l)sinh(IV^) 

z{zsmh{^/z) + -yzcosh(-yz)) 

This is Test 2 in (Duffy 1993), where it was noted that the apparent square root sin- 
gularity can be removed by expanding the hyperbolic functions. We also consider 
Test 3 in (Duffy 1993), namely 

This transform has a pole at 2 = 0, an essential singularity z = — 1/c as well as 
branch point singularities at z = and z = — 1 . With appropriate definition of the 
branch cuts, this transform is analytic everywhere off the negative real axis. 

Numerical results for transforms F2(z) and F^{z) are presented in Figures 5 and 
6, respectively. Again we see empirical confirmation of the theoretical 0{e~^-^^^^) 
convergence rate and the effectiveness of the roundoff control strategy. A compari- 
son with the results presented in (Duffy 1993) shows a faster rate of convergence 
here, which reduces the total number of nodes in the midpoint sum (4) required 
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Figure 4. Relative errors in the 00-norm when exp(— At)uo is approxi- 
mated by applying the midpoint sum (4) to the contour integral (20). For 
smaller values of N the contour (15) was used, while for larger N the 
roundoff control procedure of Section 3 was used. The three columns cor- 
respond to three sizes of A, while the two rows correspond to two values 
of t. The dashed line represents the theoretical predicted convergence 
rate 0{e-^-^^^^). 



for, say, ten-digit accuracy at f = 1, from about 26 and 36 for F2 and F3, to 18 for 
both transforms. 

Note that F^(z) is an example for which the error analysis of Section 2 can fail. 
The reason is its rapid decay as z ^> +00, particularly if r ^ 1; recall the discussion 
below (8). Nevertheless, for < r < 1 we get error curves of similar accuracy to 
those shown in Figure 6 for all values of t considered. 



4.2. Singularities on the positive real axis. When the singularities are located on 
the positive real axis, the frequency shift formula 



7(0 



-^ / e"'F(z + s)dz, 
Ini Jc 



(21) 



can be considered. The real shift s should be sufficiently large such that all the 
singularities are on the negative real axis. One may lose accuracy if s and t are 
large because roundoff errors are amplified by the factor e"^^ but in our experiments 
this did not appear to be a major factor In addition, for sufficiently large values of 
st, overflow can occur in the evaluation of the exponential factor 
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Figure 5. Relative error in the inversion of transform F2{z). For smaller 
values of N the contour (15) was used, while for larger N the roundoff 
control procedure of Section 3 was implemented. 



To demonstrate this shifting formula, we examine the transform fi(z; —5) as 
defined in (19) with pole at z = 5. Applying the method to the unshifted transform 
fails, but by shifting the singularity to we obtain the expected high rate of 
convergence, as shown in Figure 7. 

4.3. Singularities off the real axis. The theoretical predicted convergence rate of 
Q^g-i.358N^ is often observed even for moderate values of N when the transform 
has singularities off the real axis. For these calculations, however, the roundoff 
error control discussed in Section 3 does not apply and we have omitted it in the 
results presented in this section. 

The first transform of this type is Test 4 given in (Duffy 1993), 

¥4(2) = -exp f -2arcosh(\/l+z2 + z4/16) j, 

which we have simplified to 

Hz) = — \ ^. 

z(Vz2 + zVl6 + Vl + z2 + z4/16)2 

There is a pole at z = and six branch points on the imaginary axis, the outermost 
two of which are z = ±4f. The results for this transform are shown in Figure 8. 
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Figure 6. Same as Figure 5 but the transform is F^{z), with c = 0.4 and 
r = 0.5. 



We also consider Test 5 given in (Duffy 1993), namely 



2 — \/z^ — c^ 
zs/i?- — c^\J z — X^/z^ — c^ 



H^) = 



with c = (1 — X)/X and X < 1. This transform has branch point singularities on 
the positive real axis at c as well as on the imaginary axis at ztcX/ \/X^ — 1. To 
deal with the singularity on the positive real axis we shifted the transform by a 
distance c. If X = j, then c = 1 and the singularities are at 1 and ±z/\/3 before 
the shift, and at and — 1 ± f 7^3 thereafter. With proper attention to the branch 
cuts in the coding, this leads to the results shown in Figure 9. 

Despite both transforms F4(z) and F5(z) having singularities off the real axis, 
we see that the modified Talbot method still shows a high order of convergence, 
asymptotically similar to the theoretical 0{e^^'^^^^). The implied constant is large, 
however, and increases with t. In the case of Fi{z) the contour used here was not 
successful for large t. For this example the contour parameters of (Duffy 1993, Murli 
and Rizzardi 1990) are better, because they adjust for singularities off the negative 
real axis. 

For F^{z) we observe good convergence in Figure 9, except for the roundoff 
error issue. Note that in (Duffy 1993) it is mentioned that Talbot's method "does 
rather poorly" for this example. This might be due to not shifting the transform 
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Figure 7. Same as Figure 5 but the transform is Fi(z;A) with A = —5. 
The frequency shift formula (21) is used, with s = —A. 
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N N 

Figure 8. Relative error in the inversion of transform Fi{z). The contour 
(15) was used for all N shown (i.e., no roundoff control). 



or perhaps a branch-cut transgression. The results shown Figure 9 reach ten-digit 
precision with 24 nodes (i.e., 12 transform evaluations). 

Finally, we remark that in the case of singularities in the right half -plane off the 
real axis, the frequency shift formula (21) can also be used to shift singularities to 
the left half-plane. 
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Figure 9. Same as Figure 8, but for the transform F5(z), c = 1. The 
frequency shift formula (21) was used, with s = c. 



5. Software 

We have created a MATLAB function, Modif iedTalbot, which can be down- 
loaded from the web site (Dingfelder and Weideman 2013). It implements the 
method (4) on the contour (15). In its basic form the function takes as input the 
scalar transform f (z) and a single value of t. It outputs an approximation to f{t), 
and an estimated relative error. In addition, the user can 

• specify a relative accuracy (ten-digit accuracy is the default), 

• increase the maximum value of N (100 is the default), 

• specify whether the roundoff error control strategy of Section 3 should be 
employed (recommended only when the singularities of the transform are 
on the negative real axis), 

• include a frequency shift (based on formula (21)), 

• turn on or off a graphical output of the convergence history, and 

• specify an exact value for f{t) (if known), for testing purposes. 

More details on these options are given in the help lines of Modif iedTalbot. 

The algorithm starts with a small value of N, and then increases its value in 
steps of 2 until an estimated error drops below the specified accuracy. The default 
stopping test is 



Mt)-fN-2{t) 



^ 10 



-10 
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where /n(0 represents the approximation computed with N nodes in the midpoint 
sum (4). The left hand side is a reasonable estimate for the actual relative error 
when the convergence is exponential, as in f{t) — /n(0 = 0{e^ ). In the ideal 
case (singularities on the negative real axis, F(z) not exponentially small) the 
largest value of N required for this accuracy is t5rpically around 20. 

The strategy of increasing N until the desired accuracy is reached is somewhat 
inefficient, as transform values computed for /n-2(0 cannot be re-used in the 
computation of /nC^ )■ ^^ ^{^) i^ ^^t too expensive to evaluate this is not a major 
drawback, however. Relatedly, keep in mind that this code does not require any 
information on the singularities, unless the frequency shift formula (21) is used. By 
contrast, software like (Murli and Rizzardi 1990) requires the user to specify the 
locations and type of some of the singularities of the transform. We consider the 
additional transform evaluations a relatively small price to pay for the convenience 
of not having to deal with a singularity analysis of the given transform. 

The following MATLAB code snippet shows the usage of Modif iedTalbot. We 
chose test F3(z), in order to demonstrate the correct coding so that the branch-cut is 
located on the negative real axis. Also note the dot-multiplication and -division in 
the definition of the transform because the code in Modif iedTalbot is vectorized. 

» t = 1; F = Q(z) (l./z) .*exp(-0.5*sqrt(z) .*sqrt(l+z) ./sqrt(l+0.4*z)) ; 

» [f ,err] = Modif iedTalbot (F,t) 

» f = 0.722835907109632 

» err = 5. 043215570091767e-012 

For this example the exact value is /(I) = 0.7228359071, rounded to ten digits, so 
the default accuracy has been achieved. The actual relative error is about 1.8 x 
10^^'^, so Modif iedTalbot overestimates it by less than an order of magnitude. 

In Table 1 we present a comparison of Modif iedTalbot with the results from 
(Murli and Rizzardi 1990). The table shows the total number of nodes, N, in the 
quadrature sum (4), to reach an error of 10^^. In (Murli and Rizzardi 1990) this 
error was defined as the relative error if the reference value is greater than one, and 
absolute error otherwise. Note that the numbers in parentheses shown in Table 1 
are double the numbers shown in Tables 1 and 11 of (Murli and Rizzardi 1990), 
because of the N vs. 2N convention mentioned below (4). 

The first three transforms in Table 1 have singularities on the negative real 
axis. For each of these transforms the value of N required by Modif iedTalbot 
is significantly less than the number required by the algorithm of (Murli and 
Rizzardi 1990); in one instance the value of N is more than halved. The last three 
transforms in Table 1 have singularities in the complex plane. The last one of these 
has a singularity in the right half -plane, which we have shifted to the imaginary axis 
by choosing s = 1 in (21). For these three cases Modif iedTalbot yields reasonable 
results for small values of t. As t increases, however, the contour (15) moves closer 
to these singularities and the number of nodes has to be increased accordingly. 
For the cases marked (a) in Table 1 no convergence to the desired accuracy was 
achieved for N less than the default value of 100. In many of these cases, no 
convergence could be obtained even when this default value was increased. In 
these situations the algorithm of (Murli and Rizzardi 1990) is superior, but it too 
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needs large numbers of transform evaluations to reach the accuracy goal. For these 
cases it is probably wise to consider alternatives to Talbot's method. 

Table i. Comparison of the code Modif iedTalbot with the numerical 
results of (Murli and Rizzardi 1990). The numbers without (resp. inside) 
the parentheses are the values of N used by Modif iedTalbot (resp. the 
algorithm of (Murli and Rizzardi 1990)) to reach a specified error of 10^^. 
Cases where Modif iedTalbot could not reach this accuracy for N less 
than the default value of 100 are marked by (a). 



Transform 


t = 10-1 


f = 100 


t = ioi 


t = 10^ 


t = 103 


_2 

Z 


18 (22) 


18 (22) 


16 (22) 


16 (22) 


16 (22) 


log(z)/z 


14 (22) 


18 (22) 


14 (22) 


14 (22) 


14 (22) 


exp(-4yz) 


16 (22) 


14 (22) 


12 (22) 


12 (22) 


10 (22) 


arctan(l/z) 


16 (22) 


20 (22) 


46 (38) 


(a) (140) 


(a) (1050) 


log(i^) 


18 (22) 


22 (22) 


78 (46) 


(a) (224) 


(a) (5284) 


z2(z3 + 8)-i 


18 (22) 


26 (22) 


80 (46) 


(a) (304) 


(a) (2224) 



Table 2. Performance of Modif iedTalbot on six other transforms. The 
value of N to reach a specified relative accuracy of 10-^" is tabulated. Un- 
successful cases are denoted by (a) or (b), as discussed in the penultimate 
paragraph of Section 5. 



Transform 


Parameters 


t = 10-2 


t = 10-1 


f = 1 


f = ioi 


t = 102 


exp(-fcv'2)/z 


k = l 


40 


24 


22 


20 


20 




k = 5 


(b) 


(b) 


26 


22 


20 


exp(-fc/z)/z 


k = l 


20 


22 


24 


28 


44 




k = 5 


22 


22 


28 


38 


70 


y^/(z-fl2) 


a = 1 


20 


20 


20 


20 


20 




fl = 5 


20 


20 


20 


20 


(c) 


{^z + aVz + b)-' 


fl = -3, & = 4 


20 


18 


20 


20 


20 


Vz + a-y'z + b 


fl = -5, fe = 1 


18 


20 


22 


22 


24 




a = 1 


20 


22 


28 


64 


(a) 


(x/z2 + fl2)-i 




a = l 


20 


22 


34 


(a) 


(a) 




fl = 10 


22 


28 


64 


(a) 


(a) 



For additional performance tests, this time with a smaller error tolerance IQ-^*^, 
we offer six more transforms in Table 2. The inverses can be looked up in tables 
such as (Abramowitz and Stegun 1970, Sect. 29). The first five of these transforms 
have a combination of poles, branch-cuts or essential singularities on the real axis. 
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while the last one has a branch-cut on the imaginary axis. In those cases where 
the singularities are located on the positive real axis, a shift equal to the largest 
real singularity was implemented. 

In most of these cases Modif iedTalbot performs reliably, providing the default 
ten-digit accuracy with around 18 to 24 nodes in the quadrature sum (4). The 
cases of non-convergence for N less than the default value 100 are marked either 

(a) the singularities are too far off the negative real axis and t is too large, or 

(b) the transform decays too rapidly for z in the right half-plane and t is too 
small. (Note that (a) was already encountered in Table 1.) Whenever one of these 
exceptional situations occur, the graphical display of the convergence history 
offered by Modif iedTalbot is useful for identifying that something is wrong and 
that other methods should be used instead. 

Considering the results of Tables 1 and 2, we conclude that in the ideal case 
{F{z) with singularities on the negative real-axis and not exponentially small) the 
code Modif iedTalbot performs reliably. In other cases it can still perform well, 
unless one encounters one or more of the exceptions (a) and (b) listed above. We 
emphasize that the main advantages of Modif iedTalbot are rapid convergence 
and the fact that minimal singularity information has to be specified by the user. 
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